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ABSTRACT 



Ribas and collaborators have recently proposed that an additional, ~5 M e planet orbits the transiting planet host star GJ436. Long- 
term dynamical interactions between the two planets leading to eccentricity excitation might provide an explanation for the transiting 
planet's unexpectedly large orbital eccentricity. In this paper we examine whether the existence of such a second planet is supported 
by the available observational data when the short-term interactions that would result from its presence are accounted for. We find that 
the model for the system suggested by Ribas and collaborators lead to predictions that are strongly inconsistent with the measured 
host star radial velocities, transiting planet primary and secondary eclipse times, and transiting planet orbital inclinations. A search for 
an alternative two planet model that is consistent with the data yields a number of plausible solutions, although no single one stands 
out as particularly unique by giving a significantly better fit to the data than the nominal single planet model. We conclude from this 
study that Ribas and collaborator's general hypothesis of an additional short-period planet in the GJ 436 system is still plausible, but 
that there is not sufficient evidence to support their claim of a planet detection. 

Key words, stars: individual: GJ 436 - planetary systems 



1. Introduction 

The GJ436 system is unique among the nearly 250 extrasolar 
planetary systems identified so faiQ. It contains a "Hot Neptune" 
planet (planet "b") that was origi nally discovered with high pre- 
cision Doppler spectroscopy by Butler et al. (2004) and that was 
later found to transit by iGillon et al.l (l2007bl) . This planet is 
the only known member of its class that transits its host star. 
Therefore, it is an interesting target for the particular investi- 
gations that are feasible for transiting exoplanets (for a recent 
review of the observational techniques applicable to transiting 
exoplanets see Charbonneau et al. 2007). Since the discovery of 
the planet's transiting nature, follow-up studies have been car- 
ried out with the Spitzer S pace Telescope (herea f ter referred to 
as Spi t zer for brevity ) by [Gillon et ail d2007ah. iDeming et al.l 



(2007b. iDemorv et all (120071) . and ISoufhworthl (120081) : and the 
Hubble Space Telescope (HST) bv lBean et aiT(l2008l) . 

The GJ 436 system is also unusual because planet b has a 
significantly non-circular orbit despite its proximity to its host 
star. The planet has an orbital eccentricity e = 0.15 + 0.01 as de- 
termined from analyzing the radial velocities of the host star 
with the constraint provided by the observed time of the planet's 
secondary eclipse that wa s observed with Spitzer (Demin g et al.l 
120071: IDemorv et aT1l2007l) . An elliptical orbit for such a short- 
period planet (P = 2.6d) is potentially at odds with the predic- 
tions of tidal theory, which suggests that the planet' s orbit should 
become circularized on a timescale of < 10 8 yr dManess et al.l 
120071: iDeming et alj2007l) . 

Maness et al. | d2007l) found a significant linear trend in the 
radial velocities measured for GJ 436 over 6 years superimposed 



1 A regularly updated list of reported exoplanets can be found at 
http://exoplanet.eu 



on the signal from planet b. This discovery was interpreted to 
mean that GJ436 likely has an additional, but n ot necessarily 
lanet ary-mass, companion in a long-period orbit. iManess et al.l 
2007) investigated whether a long -period planet consistent with 
the radial velocity trend could be the perturber leading to exci- 
tation of planet b's orbital eccentricity. They found that it was 
possible, but far from certain owing to the unconstrained nature 
of t he object causing the slow acceleration of GJ 436. For exam- 
ple, IManess et al.l (|2007) suggested that a roughly Saturn mass 
planet in a 25 yr orbit with e ~ 0.2 would be consistent with the 
radial velocities and provide the necessary regular perturbations. 

Recently, Rib as et all d2008l hereafter RFB) have proposed 
another explanation for GJ436b's high eccentricity. They sug- 
gest that an additional, short-period planet in the system would 
provide the necessary regular dynamical impulse to the tran- 
siting planet so that it would maintain its high orbital eccen- 
tricity in the face of tidal circularization over long t imescales 
(altho ugh this result has recently been questioned by Mardling 
2008). Such a planet also met their requirement to cause the tran- 
siting planet's orbital inclination to change by 0.1 ° yr' 1 . They 
saw th is change in inclination as the reason why iButler et all 
(2004) didn't discover that GJ436b transited despite ostensibly 
having achieved the necessary photometric precision and sam- 
pling in their search. Their hypothesis is that the planet simply 
wasn't transiti ng at that epoch, whi le it was when observed at a 
later epoch by Gill on et all (l2007bl) and the subsequent investi- 
gators mentioned above. 

With the possible existence of another short-period planet 
in mind, RFB studied the radial velocities of GJ436 provided 
bv lManess et al.l d2007l) . They identified a low-significance peak 
(20% false alarm probability) in a periodogram of the single 
planet residuals and used that as the starting point for a two 
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planet fit. Assuming Keplerian orbits, they were able to obtain 
a two planet model that fit the radial velocities significantly bet- 
ter than a single planet model. 

The second planet in the RFB model has an orbital pe- 
riod, P = 5.1859d, which puts it close to a 2:1 orbital resonance 
with the transiting planet, and minimum mass, M sin i = 4.7 M m . 
RFB proposed that this second planet exists based on the suc- 
cess of their model for providing a source for the transiting 
pl anet's eccentricity, a reason for the non-detection of transits 
by iButler et al.l (120041) . and a fit to the radial velocities. If con- 
firmed this would be the lowest mass planet yet found around a 
nearby, main sequence star. Therefore, the claim deserves further 
scrutiny. 

One particular aspect of the RFB study that merits further 
investigation is the possible sensitivity of observables to grav- 
itational interactions between the two planets in their model. 
These two planets would be in close, moderately eccentric, and 
possibly non-coplanar orbits and so their mutual perturbations 
might be significant on short timescales in addition to the long 
timescales that RFB only considered. If they are, then RFB's 
claimed detection might be spurious because their model did 
not account for them. The critical issue is that Keplerian or- 
bits, which RFB used for their modeling of the radial veloc- 
ities, are only strictly valid for the two-body problem (i.e. a 
single planet orbiting a single star and also in the absence of 
significant General Relativity effects). Such orbits are only a 
sufficient approximation for modeling the data for multi-planet 
systems when the planet-planet interaction timescales are much 
longer than the length of the observations. For systems where 
short-term interactions are occurring, or are even possible, a 
model based on direct integrations of the equations of motion 
(i.e. Newtonian orbits) should be calculated to check the valid- 
ity of the Keplerian approximation. If the observables of interest 
would be significantly different when accounting for the interac- 
tions, then the Newtonian orbit model must be used. 

In this paper we assess the consistency of RFB's model 
of the GJ 436 planetary system with the observed host star ra- 
dial velocities, transiting planet primary (transit) and secondary 
eclipse times, and transiting planet orbital inclinations when us- 
in g Newtonian rather than Kep leria n orbits. Pioneering work 
bv lLaughlin & Chambers! (1200 ll) and iRivera & Lissaueri d200ll) 
have demonstrated that radial velocities with precisions on the 
order of a few ms~' are sensitive to sh ort-term dynamic al in- 
teractions for certa i n exo planet systems. lAgol et al] d2005l) and 
iHolman & Mu rray (2005) have shown that transit timings mea- 
sured with precisions of a few seconds up to a few minutes are 
quite sensitive to additional planets with masses down to even 
the terrestrial level. Transit timings for a planet near to a low- 
order resonance, as RFB propose for GJ436b, are parti cularly 
sensitive to very low-mass planets (ISteffen & Agolll2005l) . 

2. The model 

To generate a self-consistent planetary system model for 
comparing to observational data we used the Mercury code 
( I Chambers! 1 1 9991) to integrate the equations of motion. We as- 
sumed all the bodies were point masses and the only force con- 
sidered was Newtonian gravity. We chose the Bulirsch-Stoer in- 
tegration option as tests indicated this conservative method was 
necessary to achieve the desired accuracy. 

The results of the integrations were recorded with a timestep 
of 0.01 d and we used spline interpolation in the output grid of 
data to calculate model observables at arbitrary times. The com- 
bination of the selected sampling rate and interpolation method 



Table 1. Data from GJ 436b eclipses 



Parameter 


Value 


Source 


Transit time (HJD) 


2454222.61564 ± 0.00060 


1 


Transit time (HJD) 


2454280.78167 ± 0.00011 


2 


Transit time (HJD) 


2454439.41607 ± 0.00068 


3 


Transit time (HJD) 


2454444.70385 ± 0.00093 


3 


Transit time (HJD) 


2454447.34757 ± 0.00080 


3 


Transit time (HJD) 


2454463.20994 ± 0.00089 


3 


Transit time (HJD) 


2454468.49911 ± 0.00094 


3 


Secondary eclipse time (HJD) 


2454282.33 ± 0.01 


4 


Inclination^ (°) 


86.40 ± 0.10 


1 


Inclination^ (°) 


86.38 ± 0.10 


2 


Inclination^ (°) 


86.32 ± 0.08 


3 



References: (1) Re-analysis of iGillon et all d2007bf) light curve; (2) re- 
analys is of Spitzer data; (3) re-analysis of HST data; (4) Deming et al. 
d2007h . 

" At epoch 2454222.61564 

* At epoch 2454280.78167 

c At epoch 2454455.27924 



yielded accuracies of better than 0.02 ms -1 and 0.2 s in the stel- 
lar radial velocities and planet transit times respectively for test 
cases. This level of accuracy is adequate because it is more than 
an order of magnitude better than the uncertainties in the ob- 
served data that was modeled. 



3. The data 

3.1. Radial velocities 

The radial velocitie s for GJ436 we analyzed come from 
iManess et al.l (I2007I) . These same velocities were considered 
by RFB in their analysis using Kepler ian orbits. We followed 
the suggestion of IManess et al.! ((2007) and added 1.9 ms -1 in 
quadrature with the reported errors to account for additional un- 
certainties arising from stellar and instrumental sources. Because 
GJ 436 hosts a transiting planet its radial velocities wi ll exhibit 
the Rossiter-McLaughlin (RM) effect during a transit dRossiter 
ll924tlMcLaughlinll924l:lGaudi & Winnl2007l) . Therefore, these 
data cannot be used for orbit determination without including an 

additional model for the RM effect. 

We searched the Mane ss~et al.l d2007l) data for points that 
were possibly obtain ed during transit u sing the orbital ephemeris 
for planet b given bv lBean et al.l (|2008). We found that the obser- 
vations with time stamps of 2453196.772 and 2453841. 887HJD 
were obtained during a transit and we did not include these data 
in our analysis. The final data set of radial velocities we used 
contains 57 measurements spanning 6.5 yr and having a typical 
uncertainty of 3.0 m s . 

3.2. Transiting planet parameters 

We also include in the dynamical analysis a combination of pre- 
viously published and newly determined eclipse times and incli- 
nations for GJ 436b. We focused on the data that can be extracted 



from the transit light curves 
those obtained w ith Spitzer dGillqn et al 



given by iGi on et al l d2007b|) and 



2007) and HST dBean et alj 



2007a!: iDeming et all 



2008). Analyses based on differ- 



ent techniques and with different assumptions have been carried 
out on the transit light curves individually. We re-analyzed the 
photometry collectively with the same technique in order to ob- 
tain internally consistent data for the subsequent test of the RFB 
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Table 2. Single planet model parameters for GJ 436 



4. Analysis 



Parameter 



Value 



P(d) 
T c (HJD) 
T p (HJD) 

K(m s->) 
e 

wO 
i O 

<iv/<if (ins" 1 yr -1 ) 
* 2 

Degrees of freedom 

Radial velocity rms (ms 4 ) 

Transit time rms (d) 

Secondary eclipse time difference (d) 



2.64390 ± 0.0000056 
2454280.78168 ± 0.00011 
2454280.221 ± 0.083 
18.31 ±0.57 
0.140 ±0.007 
357.6 ± 11.7 
86.36 ± 0.05 
1.30 ± 0.27 
98.6 
60 
3.8 
0.00045 
-0.008 



modefl- This step has the additional benefit of increasing the 
precision on the parameters of interest, which are time varying, 
by reducing their correlation with the physical parameters of the 
star and planet that the light curves are sensitive to, which do not 
vary with time. 

We used t he e xact analytic formulae given by 
Man del & Ago! (l2002h to model the transit light curves 
and a Ma rkov Chain Monte C arlo (MCMC) method similar 
to that of iHolman et all d2006l) to identify the best fit model 
parameters. We assumed the planet and star radii were the 
same for each transit observation, but allowed each to have 
a unique central transit time and inclination. The HST data 
include observations of five partial transits spread over 11 
orbits of the planet and we determined a transit time for each 
one, but only one inclination for the group. There were 18 free 
parameters in total. We adopted the limb darkening constants for 
each light curve suggested by the authors in the corresponding 
papers initially presenting the data. We held fixed the transiting 
planet's orbital period to that determined by iBean et al.l (2008, 
2.643902 d), and its orbital eccentricity, longitude of periastron, 
time of periast ron, and velocity se miamplitude to those values 
determined by Deming et al. (2007). We assume the mass of the 
star is 0.44 + 0.04 Mq, where the uncertainty in this value does 
not significantly contribute additional uncertainty to the transit 
times and inclinations. 

The identified parameter values and errors from the MCMC 
analysis were the median and 68% confidence intervals respec- 
tively for aggregate of the trimmed chains. The transit times and 
inclinations used for testing the RFB model are given in Table Q] 

In addition to the re-determined transit parameters, we 
adopted the s econd ary eclipse time for pla net b given by 
Dem ing et alj (120071) and Demor v" st al.l (120071) . Those authors 
analyzed the same Spitzer data but reached different conclu- 
sions regarding the uncertainty in the determi ned secondary 
eclips e time. We adopted the uncertainty given bv lDeming et al.l 
(200 7D, which is an o rder of magnitude larger than that given by 
Dem orv et al.1 (120071) . The secondary eclipse times reported by 
both groups are consistent within this larger range. The value we 
used in our analysis is given in Table[TJfor completeness. 



2 We chos e to re-analyze the version of the Spitzer reduced data 
presented by iGillon et alj J2007al) with the corrected time stamps (M. 
Gillon private communication, 2008). 



4.1. Evaluating the proposed Super-Earth planet 

Our evaluation of the RFB two planet model for the GJ436 
system is based on the assumption that including an additional 
planet in a model for a system should yield a significantly better 
fit to the observational data for that system. To create a base- 
line for this test we first fit the data with a single planet model. 
We used a MCMC technique to determine the transiting planet's 
orbital parameters, a radial velocity linear trend, and a radial ve- 
locity offset to account for the relative nature of the radial veloc- 
ities. There were eight free parameters in total. The parameters 
giving the best fit model and lcr confidence intervals from the 
resulting parameter distributions are given in Table [2] The^ 2 of 
the best fit was 98.6 for 60 degrees of freedom. The larger than 
expected value for the corresponding reduced^ 2 (1.6) is due to 
some small inconsistencies between the fit and the radial veloc- 
ities, which is not unusual. The model predicted eclipse times 
and inclinations all fall within the uncertainties of the observed 
data. 

With the fit-quality baseline established we then calculated 
the dynamical model radial velocities, eclipse times, and inclina- 
tions for the orbital parameters suggested by RFB and compared 
them to the observational data. We assumed that the parameters 
were osculating orbital parameters for the mean epoch of the 
radial velocities (HJD = 2453002.7). In addition to Keplerian 
orbital parameters, RFB solved for time varying components to 
the transiting planet's period, eccentricity and longitude of pe- 
riastron. We did not fix these parameters in the model because 
if they are physical, then the dynamical calculations will natu- 
rally reproduce them. For a first test we just calculated the dy- 
namical model with the nominal parameters suggested by RFB 
and determined the radial velocity offset that gave the best fit to 
the observations. We assumed that the inclination of the transit- 
ing planet was the average of the observed inclinations and the 
second planet was coplanar. The resulting x 1 was 4.4e6, which 
is five orders of magnitude larger than the reference one planet 
model. The rms of the radial velocity and transit time residuals 
was 11.5ms _I and 0.3 d respectively. 

The orbital parameters given by RFB are not exact so we 
searched over the parameter space bounded by their given un- 
certainties while allowing the known planet's orbital parameters, 
second planet's eccentricity, velocity trend, and velocity data off- 
set to vary. Using a combination of grid search and local mini- 
mization techniques we identified a best fit model that gave a 
X 1 of 144.9. This model had the second planet non-coplanar by 
67°, which is unlikely to be a stable configuration. The smallest 
found x 1 f° r a coplanar model was 1.5e6. From this study we 
conclude that the specific two planet model proposed by RFB is 
completely inconsistent with the observational data. 



4.2. Limits to additional planets 

While the currently available data can be used to rule out RFB's 
specifically proposed Super-Earth planet at very high confi- 
dence, what do they suggest regarding RFB's general hypoth- 
esis? To answer this question we examined how well the data 
could be fit by a two planet model where the mass and semima- 
jor axis for the second planet was in the range RFB suggested 
was necessary to provide sufficient dynamical perturbations to 
planet b so that it maintains its orbital eccentricity in the face 
of tidal circularization over long timescales. We divided the re- 
gion of interest in the mass - semimajor axis plane of the sec- 
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a (AU) 



Fig. 1. Intensity map of the x 2 surface for a two planet model 
fitted to the observed data. The shading indicates the difference 
in^ 2 from the best fit single planet model. The dashed line gives 
the 3<x upper limits for the mass of a second planet as a function 
of semimajor axis. The solid line delineates the lower mass limit 
for a second planet that RFB calculate can provide sufficient per- 
turbations to excite the eccentricity of planet b. The circle shows 
the position of RFB's purported discovery. 



ond planet into 57 sub-regions. For each of these sub-regions 
we ran a local minimization algorithm initiated with some ran- 
domly selected orbital parameters for the second planet 100 dif- 
ferent times and collected the results. For each run the orbital pa- 
rameters for planet b were initiated at their single planet model 
best fit values given in Table [2] The mass and semimajor axis 
of the second planet were selected from, and restricted to, the 
sub-region limits. The eccentricity of the second planet was se- 
lected from the range 0.0 - 0.3. The argument of periastron and 
mean anomaly of the second planet were selected from their 
fully defined ranges (both 0° - 360°). The second planet's in- 
clination was selected out of, and restricted to, the range +15° 
from planet b's inclination. The smallest^ 2 found from the 100 
different runs was taken to be representative of the best fit in the 
corresponding sub-region. We carried out 500 different runs in 
three of the sub-regions to verify that only 100 iterations were 
sufficient to reliably locate the best fit. The resulting fit-quality 
map of the second planet mass - semimajor axis parameter space 
is displayed in Figure Q] 

The results from this investigation indicate that there is a 
large range of possible parameters for a hypothetical second 
planet that give a better fit to the data than the single planet 
model. The identified best fit model has x 1 - 78.7 for 54 de- 
grees of freedom. The second planet in this model has a mass 
M c = 5.0 M e , and semimajor axis a = 0.043 AU. The false alarm 
probability (FAP) for the fit-quality improvement in this case is 
5%. 

Examined in isolation the best fit solution could be consid- 
ered as evidence to support RFB's specific claim for the dis- 
covery of a second planet in the GJ 436 system because of its 
similarities with their model and in spite of the differences in the 
other orbital parameters. However, the consideration of planet- 



planet dynamical interactions has reduced the significance of the 
fit-quality improvement from that originally seen by RFB when 
using a model that does not include the interactions. The large 
FAP probability we found for the fit-quality improvement when 
incorporating a second planet in the model of system indicates 
that such an addition is not warranted. Furthermore, in the con- 
text of the larger parameter range investigation, we find that the 
best fit solution is not unique. There are solutions in other re- 
gions of the parameter space that are better than the single planet 
model and are statistically indistinguishable from the best fit. 
Therefore, we conclude that there is currently insufficient evi- 
dence to say that GJ436 definitely hosts a second planet in a 
close, exterior orbit to the already known planet, although the 
existence of such a planet cannot be ruled out. 

Upper limits to an additional planet over the range of orbital 
semimajor axes considered can be set from the results of the 
fit-quality mapping by finding the point at which increasing the 
mass increases the^- 2 above that of the single planet model more 
than a certain amount. The 3<x confidence limit corresponds to 
an increase of 9 from the baseline model. This mass limit for 
a given semimajor axis is indicated in Figure Q] We find that 
the data can only rule out a second planet with semimajor axis 
similar to what RFB propose and M c <; 8 M e at high confidence. 

5. Summary 

We have shown that a Super-Earth planet like that one proposed 
by RFB could still exist in the GJ436 system, a lthough their 
specif ic claim of a detection is erroneous. Recentlv. lAlonso et alj 
(2008) carried out another investigation into the plausibility 
of RFB's proposed planet using a novel constraint on the 
allowable inclination change of the transiting planet and reached 
a conclusion similar to ours. Ultimately, more observational 
data are needed to further constrain the architecture of the 
GJ436 system and its evolutionary history. RFB's proposed 
planet would have been the lowest-mass one yet discovered 
around a nearby star so the results of this work reemphasize 
the importance of considering planet-planet interactions when 
interpreting observations of multi-planet systems. Keplerian 
orbits are still a useful approximation for modeling many such 
systems, but their appropriateness for a particular case should 
always be tested. 
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